rm(list=ls())
set.seed(12435)
library(stargazer)

load("elections_provinces.RData")

##
## Analysis
##
elect_data_full <- elect_data
elect_data70 <- subset(elect_data, distance2<7)


####
#### Approach 1: Weimar states and Prussian provinces based on centroids
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)

###
### TABLE DA.7
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_weimar_centroid <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_weimar_centroid <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])



####
#### Approach 2: contemporary states in first stage AND Weimar states and Prussian provinces based on centroids
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)

###
### TABLE DA.8
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_weimar_contemp <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_weimar_contemp <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])



####
#### Plots for FIGURE SM5.1
####

#### AFD full sample
coeffs <- c(-.081, -.058, cfs_weimar_centroid[1], cfs_weimar_contemp[1])
sts <- c(.016, .015, ses_weimar_centroid[1], ses_weimar_contemp[1])
pdf('elections_provinces_plot1.pdf', width=6.25, height=4.2)
par(mar=c(2.5,9.2,1,2))
plot(seq(-.15,.05, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
axis(side= 2,at=rev(pos),labels=c('Baseline', 'Current state FEs\nin 1st stage',
                                  'Weimar state FEs\n+ Prussian province FEs',
                                  'Current state FEs\n+ Weimar + Prussian FEs'),las=2,tick=F, cex.axis=.8)
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD 70km
coeffs <- c(-.159, -.206, cfs_weimar_centroid[2], cfs_weimar_contemp[2])
sts <- c(.062, .058, ses_weimar_centroid[2], ses_weimar_contemp[2])
pdf('elections_provinces_plot2.pdf', width=4.3, height=4.2)
par(mar=c(2.5,2,1,2))
plot(seq(-.35,.1, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD + NPD full sample
coeffs <- c(-.092, -.067, cfs_weimar_centroid[3], cfs_weimar_contemp[3])
sts <- c(.017, .016, ses_weimar_centroid[3], ses_weimar_contemp[3])
pdf('elections_provinces_plot3.pdf', width=6.25, height=4.2)
par(mar=c(2.5,9.2,1,2))
plot(seq(-.15,.05, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
axis(side= 2,at=rev(pos),labels=c('Baseline', 'Current state FEs\nin 1st stage',
                                  'Weimar state FEs\n+ Prussian province FEs',
                                  'Current state FEs\n+ Weimar + Prussian FEs'),las=2,tick=F, cex.axis=.8)
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD + NPD 70km
coeffs <- c(-.171, -.220, cfs_weimar_centroid[4], cfs_weimar_contemp[4])
sts <- c(.066, .062, ses_weimar_centroid[4], ses_weimar_contemp[4])
pdf('elections_provinces_plot4.pdf', width=4.3, height=4.2)
par(mar=c(2.5,2,1,2))
plot(seq(-.35,.1, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()
